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We study theoretically the effects of pump-spot size and location on photon condensates. By 
exploring the inhomogeneous molecular excitation fraction, we make clear the relation between 
spatial equilibration, gain clamping and thermalization in a photon condensate. This provides a 
simple understanding of several recent experimental results. We find that as thermalization breaks 
down, gain clamping is imperfect, leading to “transverse spatial hole burning” and multimode 
condensation. This opens the possibility of engineering the gain profile to control the condensate 
structure. 


I. INTRODUCTION 

The laser has long served as a prototype for phase tran¬ 
sitions in driven dissipative systems [1,|2[. While for a 
single-mode cavity the transition is mean-field-like, in a 
multimode cavity |3] spatial fluctuations are possible, en¬ 
abling non-trivial critical behavior. In the last decade, 
there has been much interest in other examples of phase 
transitions in driven-dissipative systems. In part, this 
has been motivated by experiments on polariton conden¬ 
sation [i]-[6[- However, there are also experiments on cold 
atoms in optical cavities Hj, and proposals for experi¬ 
ments in “coupled cavity arrays” [9] 1 lj of superconduct¬ 
ing circuits [T3] or hybrid quantum systems 0- There 
are also intriguing connections between the dynamics of 
these quantum systems, and the study of similar ques¬ 
tions on the dynamics of classical “active matter” [T3 |. 
as studied in photo-excited colloidal systems [l5| . These 
systems all address common questions of how a flow of 
energy through the system affects the collective dynamics 
of the system, and the emergence of spatial structures. 

Closely related to both polariton condensates and pho¬ 
ton lasers are experiments on Bose-Einstein condensa¬ 
tion (BEC) of photons 0 in organic-dye-filled micro¬ 
cavities. Unlike polaritons, these systems have no strong 
matter-light coupling and so the normal modes are non¬ 
interacting photons. However thermalization is possi¬ 
ble 0 via the dye molecules. If a photon can be ab¬ 
sorbed and emitted many times before it leaves the cav¬ 
ity, the photon gas achieves thermal equilibrium with 
the dye. Thus, by adjusting the rates of absorption 
and emission, or cavity decay, one may interpolate be¬ 
tween an equilibrium BEC, and a strongly dissipative dye 
laser [l8[. We will refer to condensation throughout this 
paper, but we present phenomena that can be interpreted 
either as lasing or BEC. Following these experiments 
many theoretical works EMU explored topics including 
equilibration, phase coherence, and photon statistics of 
the photon BEC, and later experiments studied photon 
statistics [32 ]. 

Recently, two experiments (H, (Hj studied the spatial 
profile and dynamics of the photon BEC and their de¬ 
pendence on punrp-spot size and location, observing be¬ 
havior beyond the scope of existing models. Spatial pro¬ 
files below threshold were also previously studied in 17]. 


These works motivate this paper. Studying spatially 
varying systems moves away from the domain of sim¬ 
ple “mean-field” models of lasing or phase transitions: 
spatial modes allow for non-mean field critical behavior 
at phase transitions, and for spatial decay of coherence. 
This has been explored experimentally for polaritons in 
one- [35j and two-dimensions [36j. Considering such crit¬ 
ical behavior in extended systems, theoretical work has 
shown that features beyond the equilibrium classifica¬ 
tion 37] can arise, such as new critical exponents in three 
dimensions [Hf], the destruction of algebraic order in two- 
dimensional systems [Til , and potential novel universality 
classes in one dimension [4C|. Multimode cavity systems 
i.e. spatially extended systems — also allow for trans¬ 
verse pattern formation, as has been studied in lasers Q, 
and for polaritons [4li - 144l ]. Very recently, there has been 
an experimental realization of a system of cold atoms in 
a multimode optical cavity [45[. An important distinc¬ 
tion exists between such atomic experiments, where pho¬ 
tons couple to density or spin waves of the atoms, with 
the atom number being conserved, vs exciton-polaritons 
where photons couple to the exciton itself, creating or 
destroying excitons. Nonetheless, these systems provide 
an additional complimentary perspective on the physics 
of driven-dissipative matter-light systems. 


The aim of this paper is to introduce a model capable 
of describing how the size and shape of the pump pro¬ 
file affects the spatial profile of a photon BEC. In order 
to describe the spatial profile of a condensate, a widely 
used approach is to derive order parameter equations, i.e. 
a partial differential equation determining the time evo¬ 
lution of a field 'F(r), representing the condensate order 
parameter. The equations determining the spatial profile 
of a condensate are distinct for closed (conservative) and 
open (dissipative) systems. In a closed BEC, this equa¬ 
tion is the Gross-Pitaevskii equation [46j (GPE), which 
can be written in the form iTidtif = JK[\l/]/(5\l/*. Such 
an equation conserves an energy functional E [’!'], with 
corresponding phase evolution of the order parameter. 
In contrast, for a purely dissipative system, the time de¬ 
pendent Ginzburg Landau equation [43] (GLE) describes 
irreversible relaxation, d t 'if = — r<5E[4']/(5'k*, so that the 
final state is a state of minimum energy. A classification 
of such order parameter equations has been given by Ho- 
henberg and Halperin [37| . Combining both conserva- 
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tive and dissipative terms leads to the complex GPE or 
GLE HU, widely used for polariton condensates [4§-[II|. 
In some cases, such equations can show critical behavior 
outside the Hohenberg-Halperin classification [38]. Sim¬ 
ilar equations also arise in nonlinear optics Q, where dis¬ 
persive shifts (i.e. nonlinear dielectric functions, depend¬ 
ing on the field amplitude |\E’| 2 ) compete with dissipative 
terms describing loss and gain; for example, in a class-A 
laser, the dynamics of the gain medium can be adiabati- 
cally eliminated, leading to a complex Ginzburg-Landau 
equation of motion for the field amplitude jlj. 

In this paper we make use of a different approach, con¬ 
sidering density matrix equations of motion, rather than 
order parameter equations. This is because order param¬ 
eter equations generally only crudely model relaxation 
to a thermal state - in particular, the order parame¬ 
ter normally only describes the nracroscopically occupied 
mode(s), and neglects thermal fluctuations. For many 
examples of pattern formation in nonlinear optics this is 
entirely appropriate: no thernralization occurs, and this 
is accurately reproduced by the order parameter equa¬ 
tion. However, for photon BEC, thermalization is a key 
feature of the observed behavior, and so a complete model 
should be able to explain how this interacts with spatial 
pattern formation. Extensions of order parameter equa¬ 
tions to include energy dependent gain rates have been 
developed to address this for polaritons ( 12 , 53]. By in¬ 
cluding also noise terms phenomenologically, these can 
yield thermal distributions. In this paper we instead fol¬ 
low an approach that proceeds directly from our micro¬ 
scopic model [HI, [29( ■ We show that for weak coupling 
one can derive a tractable model combining spatial dy¬ 
namics with energy relaxation. This model describes how 
the spatial profile is determined by the competition be¬ 
tween energy relaxation and loss, and can explain the 
recent experiments [33, 34]. 

The remainder of this paper is organized as follows. 
Section El describes our model of the experiments, and 
derives a master equation for the photon modes and elec¬ 
tronic state of the molecules, eliminating the fast dy¬ 
namics of molecular vibrations. From this model, we 
derive coupled equations for the population of excited 
molecules, and the populations and coherence of photon 
modes. Using these equations, section Hill discusses the 
steady state properties of the photon cloud. In Sec. IIII Al 
we first show how, far below threshold, the occupation 
of photon modes depends both on their energy (control¬ 
ling the rate of emission and absorption for that mode), 
and also the overlap between the photon mode and the 
profile of the pump. This behavior occurs at weak pump¬ 
ing because there the excitation profile of molecules fol¬ 
lows that of the pump. The same approach allows us to 
understand how the pump profile affects the threshold 
power required for condensation, discussed in Sec. IIII Bl 
Once above threshold, the profile of excited molecules is 
significantly modified by the condensed photons, via a 
kind of transverse spatial hole burning. We discuss the 
consequences of this in Sec. IIII Cl In Section [IV] we then 


turn to study the early-time transient dynamics of a con¬ 
densate following an off-center pump. We show how the 
spatial oscillations evolve due to reabsorption of cavity 
light, leading ultimately to thermalization. Finally, sec¬ 
tion [V] provides some conclusions and outlook from our 
work. 


II. MODEL 


The photon BEC system consists of dye molecules cou¬ 
pled to photon modes in an optical cavity (see Fig. Oja). 
Each molecule has a complex optical spectra, due to 
the ro-vibrational dressing of the electronic spectrum. 
Despite this, one can nonetheless consider only two 
electronic states, the highest occupied molecular or¬ 
bital (HOMO) and lowest unoccupied molecular orbital 
(LUMO). Each of these levels is however dressed by lad¬ 
der (s) of rotational and vibrational excitations of the 
molecules. As discussed in our previous work [13[If}, one 
may adiabatically eliminate the vibrational states, lead¬ 
ing to absorption and emission rates T(±<5) for photon 
modes detuned by 6 = ui — wzpl from the Zero Phonon 
Line (ZPL) of the molecule. This results in a model in 
which the electronic state (HOMO or LUMO) of each 
molecule is explicitly represented, while the effects of the 
ro-vibrational excitations appear implicitly in the struc¬ 
ture of the rates r(±(f) discussed further below. 



FIG. 1. (Color online) (a) Cartoon of model system: 
molecules are represented by two electronic states (HOMO 
and LUMO levels), dressed by ro-vibrational excitations, (b) 
Gauss-Hermite eigenfunctions of an harmonic oscillator. 


To incorporate inhomogeneous pumping we must con¬ 
sider the overlap Um(r,:) between the transverse mode 
function of photon mode m and a molecule at r.;. We 
do not include here effects of the longitudinal mode pro¬ 
file, as we consider cases where only a single longitudinal 
mode is relevant (i.e. close enough to resonance with the 
gain medium). In these cases, except for very high or¬ 
der modes, the longitudinal mode profile does not vary 
significantly between the modes, and so any effects of 
overlap between the longitudinal mode profile and the 
excited molecules can be absorbed into a constant factor 
in the definition of emission and absorption rates. For the 
transverse modes, curvature of the cavity mirrors leads 
to an in-plane harmonic trap, so that ipm( r ) are Gauss- 
Hermite functions (see Fig. Eb)) in two dimensions and 
the corresponding frequencies are harmonically spaced, 
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oj m = u> c + ( m x + m y )e , where m combines both m x and 
wiy indices. The “cavity cutoff” oj c is set by the cavity 
length. We write the master equation describing the sys¬ 
tem as two terms, d t p = Mo[p] + The bare part 

is: 

M 0 [p] = [umaLa^ p] + '^2^C[am,p\ 

m m 

+ E^^ + >p] + E^r,p ] : a) 

i i 

where C[X,p] = 2 X^pX — [X^X,p\+. The operator 
creates a photon in mode m, and we assume all modes 
have decay rate k. The electronic state of molecule i 
is represented by Pauli operators <jf’ v ’ z . In addition to 
coupling to the cavity (see below), each molecule has a 
pumping rate T-|-(rj), and a non-cavity decay rate T^ in¬ 
corporating fluorescence into all modes other than the 
confined cavity modes. Other than the inhomogeneous 
pump, M. o[p] matches Refs. p©|29|. 

The term AIi n t[/5], describing molecule-photon inter¬ 
action, can be treated at various levels of approxima¬ 
tion, according to whether we include coherence between 
different photons modes. Including such inter-mode co¬ 
herence is numerically expensive, and is only necessary 
when significant coherence exists. The numerical cost 
arises because, if we truncate the equations to consider 
N m photon modes, the full coherence matrix scales as 
N^. As discussed later, in order to keep all significantly 
populated modes (when ksT he), we need relatively 
large values of N m . Thus, with the full equations, it is 
only feasible to simulate a few hundred picoseconds of 
time evolution, far shorter than the timescale required 
to reach the steady state. In what follows, we therefore 
first present the full equations of motion, used to study 
transient dynamics, and then introduce the “diagonal ap¬ 
proximation” , providing a more efficient approach when 
inter-nrode coherence can be neglected. 

A. Fully coherent model 

We denote the most complete form of the molecule- 
photon interaction M- mX [p\ = which takes the 

form: 

m,m' ,i 

+ Ki-Sr^lalnCF-p,a m ’cr+}} +H.c.. (2) 

The complex function A'(±<5 m ), discussed next, encodes 
the molecular absorption (emission) rate vs. the detuning 
S m = co m —uJzPL between mode m and the molecular Zero 
Phonon Line (see dashed line in Fig. [2]). 

For simple molecules, K(S) can be calculated explic¬ 
itly |29l ] . Alternatively, one may use experimentally 
measured spectra F(<5), and find K(S) by analytic con¬ 
tinuation - causality requires that K(5) is analytic 



FIG. 2. (Color online) (a) Absorption and fluorescence spec¬ 
trum of Rhodamine 6G on a (a) linear or (b) logarithmic 
scale. Points are experimental data [54j (for dye in Ethylene 
Glycol), lines show the fits r(±<5) satisfying r(<5) = r(— S)e^ s 
at room temperature. Here and throughout, we plot angular 
frequency, u> = 2nf. To indicate this we write the units as 
2ttTHz. 


in the lower half plane. As noted previously 0, HE 
|22| . thermalization of photons requires that F(A) obeys 
the Kennard-Stepanov (KS) relation j55j-[57j], T(i5) = 
T(-S)e 0s . We therefore use a function F(<5), shown in 
Fig- m that fits the experimental spectra, while satisfy¬ 
ing the KS relation. The procedure used is described in 
appendix [Al This determines F(<5) up to a prefactor. We 
denote r max = max[r(5)], and discuss below how to esti¬ 
mate r max from experimental results. The function K(5) 
is then found by standard analytic continuation of r(<5) 
to the lower half plane. 

Note that while equation © includes coherence be¬ 
tween photon modes, it neglects coherence between 
molecules. This is because dye and solvent molecules 
collide frequently causing rapid dephasing. Inter-mode 
coherence will be required to understand the dynamics, 
as studied experimentally in Ref. jUJ]. It is important 
also to note that by including inter-mode coherences, 
equation ([2]) does not make the “secular approximation”, 
which discards those bath-induced terms which are time 
dependent in the interaction picture. The secular approx¬ 
imation is often introduced as a necessary condition for 
having a completely positive master equation [HI] . How¬ 
ever several recent papers suggest the secular approxi¬ 
mation can lead to incorrect predictions [HHl - tHlI ]. As dis¬ 
cussed later, the current model falls into this class: the 
experimentally observed oscillations of the photon density 
can only occur if the molecular emission produces inter¬ 
mode coherence — an incoherent state would show no 
beating. Beyond the secular approximation, there may 
be instabilities, particularly at large r max . However, for 
large r max the Markov approximation fails & For the 
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parameters we consider, our equations are stable. 

Rather than explicitly simulating the density matrix 
p, we use the master equation above to write coupled 
equations of motion for the (Hermitian) photon correla¬ 
tion matrix [n \ m ,m' = and the coarse-grained 

excitation density, /(r) = i5(r — Within 

the semiclassical approximation Q [n] m>m / and /(r) 
obey a closed set of equations. The semiclassical ap¬ 
proximation means that we neglect correlations between 
the state of the photons and the dye molecules, so that 
expectations of products of operators can be replaced 
by products of expectations. These semiclassical equa¬ 
tions can be written in a particularly compact form by 
defining a number of other quantities. We define the 
matrices [K ±\ m ,m' = S mtrn 'K(±d m ), the mode function 
matrix [SI/(r )\ m ,m' = V’m( r )Vw(r), the overlap matrix 
f = / d d rf(r)'Jf(r) and [h] m>m / = <5 m , m '(iw m - k). We 
thus write the equations: 

d t n = hn + fp 0 K_(n + 1) + (f - IjpoK^n + H.c., (3) 
dtf{ r) = -rf(n,r)/(r) +Tf (n,r)(l - /(r)), (4) 

where po is the density of molecules. Note that in the first 
equation, while n, f are Hermitian, the matrices h,K± 
are not. In the equation for the excitation density /(r), 
the total absorption and emission rates are: 

T‘ ot (n, r) = r T (r) + 25ft (Tr [$(r)nK + ] ) (5) 

rj ot (n, r) = T l + 25ft (Tr [¥(r)K_ (n + 1)]) , (6) 

incorporating also stimulated emission to and absorp¬ 
tion from the cavity modes. To simulate these equations 
numerically, we discretize /(r) on a grid of N x spatial 
points, and then use an adaptive time-step Runge-Kutta 
approach to evolve the coupled equations. 

From these equations, we are typically interested in de¬ 
riving quantities such as the photon spectrum, and the 
photon density profile /(r) which is our focus in this pa¬ 
per. These quantities can be directly extracted from the 
photon correlation matrix. The spectrum is given by the 
diagonal elements [n] mj7rl and the photon density by 

I ( r ) = ^m( r )^m'(r)[n] m , ro /. 

771,771/ 

Simulating the full state of the system would thus re¬ 
quire solving Nf n +N x coupled differential equations. One 
can reduce this requirement by noting that elements of 
[n \m,m' that are far from the diagonal are however very 
small. As discussed earlier, the value of N m required 
is however relatively large; and even with reduction to 
terms with small | m — ml |, we find it is only feasible to 
simulate short-time transient behavior, and only in one 
spatial dimension. The results presented later involved 
150ps of simulated time requiring four hours of computer 
time, while the timescale to reach steady state is of the 
order of microseconds. 


B. Diagonal approximation 

As noted at the end of the last section, the full equa¬ 
tions are too computationally costly to allow numerical 
exploration of how the steady state profile depends on 
control parameters. To overcome this limitation, we in¬ 
troduce here the “diagonal approximation”, which is ac¬ 
curate as long as coherence between photon modes is 
small, and allows numerical exploration of the steady 
state. This corresponds to using the molecule-photon 
interaction A4 int [p] = A4^ s [p] with 

771,7 

( r (< 5 m )C[a m a?, p\ + T(-<5 m )£[a^dr, p \ j (7) 

where r(±5) = 25ft[-fT(±<5)] are the absorption(emission) 
rates, and H\ is a Lamb shift from the imaginary part 
of K(±5). This Lamb shift will however be irrelevant for 
the equations of motion as discussed next. 

In this approximation, we may write closed equations 
for the populations of the photon modes, neglecting co¬ 
herences, and thus have only N m + N x equations. De¬ 
noting the diagonal elements of the correlation matrix 
as n m = [n] TOjTn , and the diagonal overlap elements as 
fm = [f ]m,m = J d d rf(r) \ip m (r) I 2 , these coupled equa- 
tions take the form: 

= PoT( fm (j^m H - 1) 

[/“v + poT((5 m )(l fm)\ Tlmi (8) 

9tf(r) = -rf({n m },r)/(r)+rf({n ro },r)(l - /(r)). 

(9) 

Note here that the second equation, Eq. m, is identical 
to that seen previously, however the total molecular ex¬ 
citation and decay rates are now written in terms of the 
diagonal populations are: 

rrcwi.r) = rj. + ^2\ipm(r)\ 2 T(-Sm){nm + 1), 

771 

(10) 

r| ot ({n m },r) = r t (r) + ^ \ip m (r)\ 2 T(S rn )n m . (11) 

771 

In all the equations above we have kept the spatial 
dimension general, writing generic wavefunctions 
The experiments are in two dimensions, in which case the 
mode functions should take the form: 

, ,, {it) H ~. (it) e ' r ’ ,2 '“ 

U>m( r) = --- - -, 

feo \J ir2 mx + m y m x \m v ! 

where we take to = (to x , m y ) as a combined index and 
H m (x ) is the mth Hermite polynomial. In certain cases, 
it is possible to efficiently find the steady states in two 




5 


dimensions, and where this is possible, we follow this 
approach. However, when directly solving the equations 
of motion, it is intractable to keep all two dimensional 
modes with energies Tiu < kpT, and so in some cases 
below we instead restrict to one dimension, for which: 




1 


Hr, 



e ~x 2 /2t 


2 

HO . 


In the following we will present analytic results for gen¬ 
eral dimension d , and specify d = 1 or d = 2 for the 
numerical results. 


III. STEADY STATE 


In the following we explore the consequences of a finite 
size Gaussian pump spot, 


r t (r) 


rf r 

(27TO-2 )d/2 eX P [ 2<Tp J ! 


where ri nt is the integrated intensity, ap the spot size, 
rp the offset, and d the dimension. Note that IT 11 * has 
dimensions of [T] -1 [L] d . Similarly, since r(±<5) are mul¬ 
tiplied by po or |7/> m (r)| 2 , this means r max also has di¬ 
mensions \T\~ 1 [L] d . We measure all lengths in units of 
the oscillator length £ho of the harmonic trap potential, 
and measure all (d-dimensional) densities in units of 
For comparison to Ref. [33| , in the first part of this paper 
we set rp = 0 . 


A. Far below threshold 


Far below threshold, when T^-(r = 0) <C T^, both n„ 
and /(r) are small, so the steady state of Eq. (18141) is 


/( r ) — 


Ff(r) 


r) ot ({n m = 0},r)’ 


— frt 


r (~6 m ) 


‘r(5 m ) + n/po' 
( 12 ) 

Note that in the denominator of the expression for /(r) 
we have not written Tj,, but rather r^ ot ({n m = 0},r), 
which includes also the spontaneous emission into empty 
cavity modes. However, for relevant parameters (see be¬ 
low), the cavity mode contribution to r) ot ({n m = 0},r) 
is small, so /(r) ~ r-|-(r)/Fj., and the overlaps f m de¬ 
pend on the pump shape. In this limit the shape of /(r) 
depends only on the shape of the pump, the normalized 
spectrum r(±<5)/r max , and the dimensionless parameter 
r) = fv/pobmax- If p -C 1 and the KS relation is obeyed 
then: n m = / m e _/3,5m . If one also has ap » £no, then 
f m is independent of m, and so there is a thermal photon 
distribution leading to a thermal photon cloud profile: 


/(r)oc^e /3<5m |'0m(r)| 2 oc exp ( - 


r 2 \ 
2 ( 7 ^ J 



0 5 10 15 

ap/t H o 


FIG. 3. (Color online) (a) Photon cloud size and (b) pho¬ 
ton number (per unit power) vs pump-spot size, for various 
77 = /t/por max . Note that in panel (a), the lines for the two 
smallest values of 77 lie on top of each other. Plotted for d = 2, 
far below threshold using the closed-form solution in Eq. m 
with u} c = 3200THz. Dashed lines in (a) show thermal size 
or (see text) and pump size ap. The points marked by the 
symbols correspond to the points for which cross-sections are 
shown in Fig. [4] 


with ap = Zho /\/2 tanh(/3e/2), which recovers the clas¬ 
sical thermal cloud size if kpT » e, see Fig. EDa). 

Thermalization fails for small ap or large 77 . This can 
be seen by looking at the actual cloud profiles, as shown 
in Fig. [4] At small ap this failure is due to the mode 
dependence of f m : A small pump spot populates only 
the low order photon modes, leading to an unnaturally 
small (i.e. cold) photon cloud, even when 77 <C 1 (see 
Fig.[4ja)). Note that for this non-thermal distribution to 
occur in this limit, the presence of the non-cavity decay 
rate Tj. is crucial: if both sources of loss and k are 
small, repeated absorption and re-emission of photons 
will occur, producing a thermal distribution. For large 77 
thermalization fails because n m ~ / m r(- S m )/Kpo, and 
so no Boltzmann factor arises. In Fig. (3])a) this gives 
a photon cloud which is larger than ap (see outermost 
(cyan) line in all panels of Fig. [4]). 

Figure EKb) shows the total photon number (per unit 
of incident power) vs pump spot size. For large spots, 
the number falls off as ( ap)~ d . This is because for all 
modes m with extent much smaller than ap, one may 
approximate 


fm 


r t (r = 0 ) 


pint 

i t 


r.|.(\/27rcrp) d ’ 


(13) 


thus giving Atotai/R,. 11 * = Yf, m fm e P Sm /Y'f 1 oc a p d . In 
contrast, for small spots, the number saturates; here ap 
is much smaller than the extent of the relevant modes 
and so f m — \4>m{ r = COpr^/r^., independent of ap. 
This expression clearly means that the occupation of all 
odd modes (which have a node at r = 0) will vanish. 
Even modes however have a value for i/j m (r = 0), and so 
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FIG. 4. (Color online) Photon cloud profile far below thresh¬ 
old for various 77 = ft/poF max . Note that in all panels, the lines 
for two smallest values of r/ are indistinguishable. The gray 
dashed line indicates the equilibrium profile, and the dash- 
dotted line indicates the profile of the pump spot. Plotted for 
op/l-ao = 1,5,9 respectively (a-c) and other parameters as 
in Fig. [3] 



FIG. 5. (Color online) (a) Photon cloud size and (b) photon 
number (per unit power) vs cavity cutoff u) c for various p = 
K /Pc>r m ax. For comparison to later figures, this is plotted for 
d=l. Spot size op = 16Ajo; all other parameters as in Fig. [3] 


iVtotai clearly saturates at a finite value. Once op <C (-no, 
this overlap with even modes becomes independent of 
op, leading to the saturation observed in Fig. [3] In ex¬ 
periment [331 ], the total photon number initially increases 
with spot size, an effect not seen here. Such a discrepancy 
could perhaps arise if the coupling of small pump spots 
into the cavity is less efficient due to some aspect of the 
pumping optics: given the above arguments about over¬ 
laps for small pumping spots, it is hard to explain such 
a low efficiency when considering purely light trapped 
inside the cavity. 

The behavior in Fig. (STa) for rj cs 10 -3 is very simi¬ 
lar to the experimental results of Ref. [331]. Using other 
known parameters of this experiment, po — 10 8 ^ho’ an< ^ 
k = 500MHz, this gives r max = 5kHz£^ 0 . We use these 
parameter values below unless otherwise stated. In or¬ 
der to verify the assumption made earlier that we may 
replace r) ot ({n m = 0},r) ~ F^, we use these values to 
estimate the effect of loss into empty cavity modes. Com¬ 
paring the rate of emission into the lowest cavity mode, 
F(-<5 m )|^ 0 (r = 0 )| 2 < r max /( v / 7 T^ 0 ) ~ 2.8kHz, to the 
observed background decay rate Tj, ~ 250MHz, this im¬ 
plies that even if the first 1000 cavity modes contributed 
to the emission equally, the cavity-mediated contribution 
would be far smaller than the background. The contribu¬ 
tion of high order cavity modes however falls off due both 
to the overlap ip m (r), and the eventual decay of T(— S m ) 
at large m. Thus, the assumptions made at the start of 
this section are indeed justified. 

So far in this section we have explored dependence on 
the properties of the pump spot. Another relatively easy 
parameter to tune is the cavity cutoff frequency w c . In¬ 
deed, as discussed extensively by [34j], tuning this param¬ 
eter can be used to control the degree of thermalization. 
A large value of ui c will enhance reabsorption of cavity 
and thus lead to thermalization, while a smaller value re¬ 
duces reabsorption and prevents equilibration. Later in 
this paper we discuss this behavior at and above thresh¬ 


old. Figure [5] shows the effect of cutoff frequency on 
properties far below threshold. One can clearly see in 
Fig. |5fa) that at large w c , the photon cloud size ap¬ 
proaches the equilibrium cloud size. Similarly, consider¬ 
ing the total number of photons, one sees that at large ui c 
the behavior asymptotically approaches the equilibrium 
result N = fme~ l3Sm indicated by the gray dashed 
line (where f m is given by Eq. (fl3l) as we consider a 
large cloud size, op ino)- At smaller to c , the photon 
number decreases as N ~ ]T) m fmX{—S m )po/K, hence the 
strong dependence seen upon the value of 77 = K/poTmax- 


B. Threshold pump power 

We next turn to the behavior at threshold, and ex¬ 
plore how this depends on the pump size. As with any 
finite size system, the threshold is not perfectly sharp; 
for definiteness we use the same threshold condition de¬ 
fined in Refs. [ 22 ], [29|. Figure [G] shows the threshold 
value of ri nt vs cavity cutoff, lo c , and vs pump spot size, 
op. These calculations, time-evolving Eqs. mm nu¬ 
merically, are computationally expensive in d = 2 so we 
consider d = 1 from hereon. In equilibrium the threshold 
condition is r-j-(r = 0) = where S c = u> c — wzpl 

(see 0). This condition has a simple meaning: it 
identifies when the effective chemical potential of the 
molecules, /r e ff(r) = iczpt, + fcBTTn[rt(r)/r.i.1 reaches the 
lowest photon mode [29(, w c . 

As has been discussed previously [22|, 0, it is 

notable that the lasing transition, normally associated 
with inversion, can be described here as corresponding to 
a thermal distribution with a positive temperature and 
chemical potential /r e ff < wzpl- The fundamental rea¬ 
son why electronic inversion is not required here is the 
different rates of absorption and emission, T(±<5). Net in¬ 
version of the electronic state is normally required for las- 
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FIG. 6. (Color online) Threshold (integrated) pump power, 
calculated in d = 1 for various values of F max . (a) vs cavity 
cutoff, ui c , (b) vs pump spot size. Simulations performed in 
one spatial dimension with N m = 200 photon modes, and 
N x = 300 spatial grid points. The dashed line shows the 
equilibrium result T^ nt = \/2TiapY^{r = 0) = \f2napT\e^ 5c . 
At large spot sizes, the threshold power increases as (crp)“ and 
so is linear for this ID simulation. Panel (c) shows the ratio of 
threshold power vs equilibrium threshold power for the region 
highlighted in panel (a), illustrating the asymptotic approach 
to equilibrium. Panel (d) shows an enlarged region from panel 
(b) as indicated. 


ing because absorption and emission coefficients match, 
so net gain requires an inverted population. Here, for 
modes with 6 = to — wzpl < 0, we have r(— S) > F(<1), 
and so for these modes, emission exceeds absorption even 
without electronic inversion. If one considers the micro¬ 
scopic ro-vibronic levels of the molecule, there is inversion 
between the lowest ro-vibrational level of the electronic 
excited state, and the higher ro-vibrational levels of the 
electronic excited states, and it is these transitions that 
have net gain [TH], However, in our model, the fast dy¬ 
namics of the ro-vibrational levels have been adiabati- 
cally eliminated. 

The thermal equilibrium prediction of threshold at 
// e ff ( r = 0) = w c is shown as the dashed line in Fig. [6] The 
actual threshold in Fig. (HKa) is however non-monotonic. 
At large w c the system is thermal, and so threshold in¬ 
creases exponentially with oj c . Figure [HKc) illustrates 
the asymptotic approach to the equilibrium behaviour by 
plotting the ratio T™*/(-\Z2napTie^ Sc ), which approaches 
1 at large oj c . At small w c the absorption and emission 
rates are too small to compete with cavity loss, and so the 
threshold pump increases. Such non-monotonic depen¬ 
dence has been seen experimentally [33| . The minimum 
of threshold becomes more pronounced as one increases 
the cavity loss rate n or decreases the peak emission rate 


r 


max* 


In Figure [61(b), we see that in d = 1, T™* oc op at 
threshold, except for small spot sizes where it saturates 
(see the enlarged region in Fig. EKd) ). From the asymp¬ 
totic form of the equations it is straightforward to see 
that in d dimensions this result becomes F^ nt oc (<Jp) d - 
Such a dependence on spot size occurs because threshold 
is reached first at the trap center, where /i e ff(r) is great¬ 
est. As such, it is the peak pump power oc F| nt tTp d that 
enters the threshold condition. It is however important 
to note that the simple power law arises only for large 
enough spot sizes, whereas for spot sizes comparable to 
the harmonic oscillator length, saturation of the critical 
integrated power occurs. In the (two dimensional) ex¬ 
periment of Ref. j33j a phenomenological power law with 
exponent ~ 1.5 was extracted from a least squares plot 
to data on a log-log scale over one decade of spot size. 
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FIG. 7. (Color online) Clamping of the gain profile for 
sufficiently large r max (in d = 1). Panels (a,b,c) show 
Fmax = 5kHz^HO, panels (d,e,f) r max = lOOkHz-feo- Top 
panels show the gain profile /(r) (solid), the value set by 
the pump, fp( r) = r^(r)/(ri + T-f(r)) (short dashed), and 
the clamped value /p = [e~^ Sc + 1] _1 for the cavity cutoff 
c d c = 3200THz. Middle panels show the (normalized) photon 
density plotted on a logarithmic scale. For comparison the 
blue dashed line shows the profile of the ground mode \i/jo(r)\ 2 . 
Bottom panels show the mode populations n m in comparison 
to an equilibrium Bose-Einstein distribution, demonstrating 
multimode condensation. Simulations performed in one spa¬ 
tial dimension with N m = 200 photon modes, and N x = 300 
spatial grid points. 




























C. Above threshold — transverse hole burning 


Once threshold is reached, in equilibrium, the chemical 
potential locks at /i e ff = w c . This means that the “gain 
profile” /(r), i.e. the fraction of excited molecules, must 
saturate, /(r) < /e = [e _/3l5c + l] -1 . Such a saturation 
of /(r) is also expected for a laser, and is known in that 
context as gain clamping. As noted above, because of the 
different rates T(±<5), no electronic inversion is required 
for lasing, and so net gain exists even though /(r) < 1. 
As such, for a photon BEC, the laser concept of gain 
clamping and the thermal concept of chemical potential 
locking are the same. Gain clamping or chemical equi¬ 
librium also imply that /(r) should become uniform at 
threshold, and we next turn to explore if and how this 
occurs. Figure [7] shows /(r) slightly above threshold for 
two values of r max . At r max = 100kHztAo> clamping is 
seen near the trap center, but for r max = 5kHztAo it is 
absent. The dependence on T max follows from the steady 
state result, 

, , Tf «n m },r) 

HT) Tf({n m },r) + r‘ ot ({n m },r) 

and the form of T|°|({n m },r) in Eq. (1101111) . If n m is a 
Bose-Einstein distribution with chemical potential /a and 
T(<5) obeys the KS relation then 

T(—<5 m )(n m + 1) = T(8 m )n m e~ 

This means that if both the following are obeyed: 

^ |'0m(r)| 2 T(—<5 m )(n m + 1) » Vi, 

m 

Y, \tp m {r)\ 2 T(5 m )n m > T t (r), 


then one has an uniform gain profile /(r) = Je- We thus 
see that gain clamping requires large T(±d). Moreover, 
since the condensed rnode(s) are concentrated at the trap 
center, gain clamping is spatially restricted, as seen in 
Fig. [ltd). In the terminology of lasers, this is analogous 
spatial hole burning [63]. In standard laser resonators, 
hole burning is discussed in terms of competition among 
different longitudinal modes, as typically laser resonators 
are much longer than the wavelength of light, but de¬ 
signed to support few transverse modes. In the photon 
BEC the situation is opposite: the microcavity supports 
only one longitudinal mode nearly resonant with the gain 
medium, but many transverse modes. As such, one has 
“transverse spatial hole burning”, leading to patterns in 
the gain medium as a function of the transverse coordi¬ 
nate, as opposed to the more standard patterns along the 
cavity axis. 

As the gain clamping is imperfect, other modes may 
reach threshold leading to multimode condensation. In 
Fig. El panels (c,f) show the mode populations n m , 
demonstrating that several modes are macroscopically 
occupied. For values of r max larger than those shown, 


multimode condensation is suppressed. Both the inho¬ 
mogeneous /(r) and multimode condensation are signa¬ 
tures of imperfect thermal equilibrium. An interesting 
question for future work is how such hole burning might 
be used to engineer the photon condensate profile. 


IV. DYNAMICS 

Having discussed the effects of the pump-profile on 
steady-state properties, we now turn to consider dynam¬ 
ics, and the transient response following a pump pulse. 
The calculations in this section are motivated particu¬ 
larly by the work of Schmitt et al. j34|, who studied the 
dynamics of the photon condensate after an off-center 
pump pulse and observed oscillations of the photon con¬ 
densate. These oscillations correspond to transverse mo¬ 
tion of a photon wavepacket in the effective harmonic 
trap potential of the mirrors. Schmitt et al. 34] showed 
in particular that for a higher frequency cavity cutoff to c 
(closer to the peak of the molecular emission and ab¬ 
sorption spectrum), thermalization occurred, while for a 
lower cutoff, oscillations persisted to later times. Schmitt 
et al. | 34 | also provided a theoretical discussion of their 
results, making use of a set of semiclassical equations 
for quantities n m ( r), i.e. spatially dependent popula¬ 
tion of a given mode. Such a quantity does not appear 
within the model discussed above: a given mode has a 
given spatial profile |'0 m (r)| 2 , and as discussed earlier, 
this means that within the diagonal model one cannot 
have asymmetric distributions, since the mode profiles 
| (r) | 2 are all even. Our aim in this section is therefore 

to show that our model can reproduce the behavior seen 
by Schmitt et al. [3T| , while considering the full covari¬ 
ance matrix [n ]m,m', or equivalently its spatial represen¬ 
tation n(r, r') = [ n ]m,m'V , m(r)V’ m (r')- 

To explore both thermalization and oscillations, it is 
crucial to use the model as presented in Eq. ©, i.e. with¬ 
out the secular approximation. This can be seen from 
quite general arguments. Firstly, in order to describe 
an off-center photon pulse, we must allow emission into 
wavepackets, not just populations of eigenstates, since 
\t/j m (r )\ 2 is symmetric for all modes, and so cross terms 
( r ) are crucial to give an off center photon dis¬ 
tribution /(r). Including cross terms is also crucial in or¬ 
der to describe oscillating wavepackets since time depen¬ 
dence occurs via beating between modes. This is incom¬ 
patible with the standard approach of secularizing mas¬ 
ter equations to produce a Lindblad form: secularization 
is appropriate if one may assume that any cross terms 
between different modes oscillate fast, and so should be 
removed. The result is an equation which can only pro¬ 
duce populations of modes. Physically it is however clear 
that the beating between modes is not a fast process to 
be eliminated, but a process on timescales comparable to 
emission and absorption. 

Naively, including cross terms between photon modes 
might suggest an alternate phenomenological equation 
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which is of Lindblad form: 


•Mint \p\ 


T{-5)C 


m,i 


+ T(5)C 


'^2ip m (r i )a m af, p 

m,i 


(14) 


However, such a form is not able to describe tliermaliza- 
tion, and its dependence on cutoff wavelength. As dis¬ 
cussed earlier, such thermalization relies on the fact that 
emission and absorption rates on the detuning of a given 
mode, to, but by its form, Eq. m has rates independent 
of mode. Modeling both the emission into wavepackets 
(i.e. inclusion of cross terms), and the mode-frequency- 
dependence of emission rates requires an equation of the 
form of Eq. ©. 
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FIG. 8. (Color online) Oscillations following an off-center 
pump pulse (in d = 1). Panels(a,b) are for uj c = 3200THz, 
where thermalization is not sufficient to suppress the oscilla¬ 
tions. Panels (c,d) are for u> c = 3250THz, showing a transi¬ 
tion to a central time-independent photon cloud at late times. 
Panels (a,c) show /(r), and (b,d) show I( r). The pump pulse 
is a Gaussian at rp = 7fao, with width up = 0.3-fao, du¬ 
ration duration 5ps and intensity T). n = 24 GHz£ho- Other 
parameters; k = 100MHz, Tj. = 250MHz, Tmax = 3 MHzHjo, 
mode spacing e = 0.4THz. These last two parameters 
are enhanced compared to experiments to reduce simulation 
time. Simulations performed in one spatial dimension with 
Nm, = 180 photon modes, keeping inter-mode coherences for 
modes with \m — m!\ < 40, and N x = 300 spatial grid points. 
This gives 15, 000 coupled equations, and this requires four 
hours to simulate 150ps. 


Using Eq. m , we simulate the dynamics following a 
short, high intensity, off-center pump pulse. The result¬ 
ing photon density I(r) and excited molecule fraction 
/(r) are shown in Fig. [5] The behavior differs according 
to the cavity cutoff oj c as seen experimentally (34|. Os¬ 
cillations occur initially in both cases, but at late times, 
they are replaced by a cloud near the trap center when 
uj c is large enough. Note that, as seen in experiment, the 


switch to the thermal condensate does not occur through 
a continuous damping of the amplitude of the oscilla¬ 
tions, but rather through a growing intensity of the cen¬ 
tral cloud, and decaying intensity of the oscillating cloud. 

As first shown by Schmitt et al. |34j, the origin of 
thermalization is that thermalization occurs when re¬ 
absorption of photons leads to a flat gain profile /(r) ~ 
fE in the center of the trap. Our model also reproduces 
this behavior, as can be seen from Fig. [Sic), and is also 
more clearly shown in Fig. [9] which plots cross sections of 
/(r) at various time slices. 



FIG. 9. (Color online) Cross sections of Fig. [8] for uj c = 
3250THz. The gray dashed in panel (a) is the equilibrium 
population of molecules. In panel (b) the blue dashed line 
shows the profile of the ground mode |i/’o(r)| 2 . 

Figure [TO] shows how the photon spectrum evolves with 
time for w c = 3250THz, where oscillations disappear at 
late times. As also seen by |s3|, the high energy modes 
rapidly reach a steady population, while the lower energy 
modes evolve more slowly, and are still evolving even af¬ 
ter 600ps of simulation time. Note however that the occu¬ 
pation of the high energy photon modes does not match 
the dye temperature. This is because of the limited spa¬ 
tial extent of the gain profile /(r) - this extends over a 
range r/t < 7, so modes up to to < 50 will be effectively 
populated. This corresponds to uj m = to c + me < 3270, 
higher modes start to be suppressed by spatial overlap, 
leading to colder photon distribution, as discussed ear¬ 
lier. 


V. DISCUSSION AND CONCLUSIONS 

In conclusion, we have presented a theoretical model 
capable of describing the spatial dynamics, relaxation, 
and thermalization of an inhomogeneously pumped pho¬ 
ton condensate. Using this, we have reproduced recent 
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FIG. 10. (Color online) Spectrum (i.e. diagonal elements of 
correlation matrix, ]n] m , m ), plotted at the same times as the 
cross sections in Fig. [9] plotted for w c = 3250THz. The black 
dash-dotted line in involves fitting a Bose-Einstein distribu¬ 
tion, the fitting temperature is 175K, and chemical potential 
/r = 3251THz — see further discussion in the text. The gray 
dashed line is for comparison the Bose-Einstein distribution 
with the Equilibrium parameters T = 300K, /i = 3250THz 


experimental results studying the effects of small pump 
spots. Even without photon loss, thermalization can be 
inhibited by small spot size. Our model gives direct ac¬ 
cess to the gain profile /(r), which is hard to access ex¬ 
perimentally. By doing so, we see the observed behavior 
at and above threshold is related to gain clamping and 
spatial hole burning. 

In this paper we have presented results only for Gaus¬ 
sian pump spots and harmonic trapping potentials. How¬ 
ever, the equations we present can easily be generalized 
to other pump profiles or trapping potentials, e.g. ring- 
shaped pumps. As such, Eqs. m or the diagonal ap¬ 
proximation, Eqs. (BED , provide a useful model to predict 
how the spatial and spectral structure of photon conden¬ 
sates is determined by the properties of the pump. As 
well as the computation framework, some general princi¬ 
ples can be identified from our results. Far below thresh¬ 
old, the condensate profile is simply given in terms of 
the overlaps f m between the pump profile and a given 
trap mode, allowing one to understand how the cloud 
will shrink if f m is significant only for low-order modes. 
A ring pump would in contrast lead to overlaps for a spe¬ 
cific range of m x ,m y , and a corresponding ring-shaped 
thermal cloud. Far above threshold, the picture of gain 
saturation and spatial transverse hole burning can pro¬ 
vide an intuitive picture of which condensate modes are 
favored or suppressed by a particular pump profile, and 
which profile shapes favor single or multimode condensa¬ 
tion. 

Mode competition and spatial pattern formation in 
driven dissipative systems have recently prompted signif¬ 
icant interest in other contexts, e.g. for random lasers §34|, 
[65j |, where the role of mode competition and the statistics 
of multimode lasing have been studied. Mode competi¬ 
tion is the basis of transverse pattern formation in non¬ 


linear optics sm, and is a prime example of pattern 
formation out of equilibrium [671 ]. The model and results 
presented in this paper provide the foundation to study 
these effects in the photon condensate. 

As noted in the introduction, for both lasers and con¬ 
densates, order parameter equations are widely used to 
describe the spatial dynamics. In contrast, our work here 
is based on solving equations for the photon density ma¬ 
trix directly. By solving all elements of the photon den¬ 
sity matrix, we allow for populations of modes in addi¬ 
tion to the condensate mode. Such a treatment is crucial 
in reproducing thermal expectations, particularly below 
threshold. In the absence of noise terms, this is not possi¬ 
ble with standard order parameter equations. However, 
thermal fluctuations can be incorporated into classical 
field methods by adding stochastic noise terms. This 
is discussed extensively in the review article by Blakie 
et al. (68]. It is an interesting question for future work to 
develop a stochastic order parameter equation that can 
reproduce the dynamics studied here. 

As compared to work on pattern formation in polari- 
ton condensates (H, jill - lil ] , an advantage of the photon 
condensate system is that we possess a clear model of the 
processes leading to the thermalization, and one which is 
readily tractable for spatially extended systems. For po- 
laritons, various phenomenological models [52j, [53} have 
been developed, and attempts to derive models micro¬ 
scopically [691 ] have been made. However questions re¬ 
main open about the relative role of polariton-polariton 
scattering vs. scattering with phonons in the semiconduc¬ 
tor fZCllZl}. In contrast, the models provided in this pa¬ 
per for the weak-coupling photon condensate are far sim¬ 
pler, as the effect of the (localized) ro-vibrational modes 
are well characterized through the function T(<5). An in¬ 
teresting question arising from this is to explore how to 
extend the treatment presented here to the case of strong 
coupling with organic molecules [73 - 174} . 
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Appendix A: Extracting F(<5) from experimental 
spectra 

As discussed in section El we aim to use the full spec¬ 
trum, P((5), extracted from experimental measurements. 
This spectrum includes the effects of all ro-vibrational 
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modes automatically. The experimental measurements 
provide two functions r abs . iexp (w), F^or., exp M, corre¬ 
sponding to absorption and fluoresence measurements. 
This appendix describes the procedure we use to find a 
spectrum consistent with both these measurements and 
the Kennard Stepanov relation. 

We first produce a single experimental function T(<5), 
by averaging the absorption and fluorescence measure¬ 
ments. This is done by identifying ujzpl from the mid¬ 
point between the peak absorption and emission, and 
then shifting and overlapping the experimental spec¬ 
trum about these points to produce an averaged exper¬ 
imental function. We use a cubic-spline fit to the ex¬ 
perimental data, and construct the function T exp ((5) = 
[rabs.,exp(^zPLT^)“t“Tfl uor; exp(^zPL <5)]/2. This yields a 
single function, but not one consistent with the Kennard- 


Stepanov relation. Furthermore, at large negative <5, 
where T exp (<5) is small, it falls below a noise floor, and 
so the experimental measurements cannot probe the ex¬ 
ponentially small values that must be present to satisfy 
Kennard-Stepanov. 

To address both the above points, we then construct 
the function: 

r(o) 2 r exp (o) + ^ F exp ( o)e , (Al) 

where x(5) interpolates smoothly from — 1 at large nega¬ 
tive S to +1 at large positive S, such that x(—5) = — x(S). 
One may readily check that this ensures T(— S)e l3S = T(<5) 
as required. The interpolation has the effect that we 
use T(<5) ~ r exp (<5) where T exp (<5) is large, and T(<5) ~ 
r exp (~ S)e l3S where T exp (<5) is small and below the noise 
floor, but T exp (— S) is large — see Fig. [2] 
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